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We estimate the strength of the bandpass-integrated thermal emission from the extrasolar planet 
TrES-4 at 3.6, 4.5, 5.8, and 8.0 ^m using the Infrared Array Camera (IRAC) on the Spitzer Space 
Telescope. We find relative eclipse depths of 0.137 ± 0.011%, 0.148 ± 0.016%, 0.261 ± 0.059%, and 
0.318 ± 0.044% in these four bandpasses, respectively. We also place a 2a upper limit of 0.37% on 
the depth of the secondary eclipse in the 16 //m IRS peak-up array. These eclipse depths reveal 
that TrES-4 has an emission spectrum similar to that of HD 209458b, which requires the presence 
of water emission bands created by an thermal inversion layer high in the atmosphere in order to 
explain the observed features. TrES-4 receives more radiation from its star than HD 209458b and 
has a correspondingly higher effective temperature, therefore the presence of a temperature inversion 
in this planet's atmosphere lends support to the idea that inversions might be correlated with the 
irradiance received by the planet. We find no evidence for any offset in the timing of the secondary 
eclipse, and place a 3cr upper limit of |ecos(a;)| < 0.0058 where e is the planet's orbital eccentricity 
and uj is the argument of pericenter. From this we conclude that tidal heating from ongoing orbital 
circulatization is unlikely to be the explanation for TrES-4's inflated radius. 

Subject headings: infrared: techniques: photometric - eclipses - stars:individual: TrES-4 - planetary 
systems 



1. INTRODUCTION 

Transiting extrasolar planets offer a unique opportu- 
nity to study the diversity of planetary atmospheres, and 
also provide an important testing ground for models of 
these atmospheres. The "hot Jupiters" are a prime ex- 
ample: although these planets have masses similar to 
the gas giant planets of the Solar system, they orbit at 
less than 0.05 AU from their parent stars. Because the 
time scale for tidal synchronization is short compared to 
the ages of these systems, these planets are expected to 
be tidally locked, with permanent day and night sides. 
This presents a significant challenge for planetary atmo- 
sphere models, as the equilibrium temperatures are sub- 
stantially higher (1000 — 2000 K) and the atmospheric 
circulation patterns significantly different than those of 
Jupiter. The hottest of these planets have temperatures 
comparable to those of the coolest stars, placing them in 
a unique parameter space with potentially exotic atmo- 
spheric chemistry. 

By measuring the decrease in flux as these planets 
move behind their parent stars, it is possible to charac- 
terize the light emitted by the day sides of these planets, 
and to construct a rough spectrum (Charbonneau et al. 
2005, 2008; Deming et al. 2005, 2006, 2007; GriUmair et 
al. 2007; Richardson et al. 2007; Harrington et al. 2007; 
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Demory et al. 2007; Knutson et al. 2008; Machalek et 
al. 2008). These observations indicate that there may be 
two distinct classes of hot Jupiter atmospheres. One class 
of planets, including HD 189733b (Deming et al. 2006; 
GriUmair et al. 2007; Charbonneau et al. 2008; Barman 
et al. 2008) and TrES-1 (Charbonneau et al. 2005), have 
emission spectra that are consistent with standard ID 
cloud-free atmosphere models for these planets (Hubcny 
et al. 2003; Sudarsky et al. 2003; Seager et al. 2005; Bar- 
man et al. 2005; Fortney et al. 2005, 2006a, b; Burrows 
et al. 2005, 2006, 2008). These model spectra are dom- 
inated by strong absorption features from CO and H2O 
in the infrared. In contrast, the emission spectra of plan- 
ets such as HD 209458b (Deming et al. 2005; Richardson 
et al. 2007; Burrows et al. 2007; Knutson et al. 2008) 
and XO-lb (Machalek et al. 2008) require models with a 
temperature inversion between 0.1 — 0.01 bars and water 
bands in emission instead of absorption (Fortney et al. 
2006a, 2008; Burrows et al. 2007, 2008) in order to ex- 
plain the observed features. The hot Neptune GJ 436b 
(Deming et al. 2007; Demory et al. 2007) and the hot Sat- 
urn HD 149026b (Harrington et al. 2007) have both been 
observed in the Spitzer 8 /im bandpass, and analysis of 
observations at additional wavelengths is pending. 

Although the nature of the high-altitude absorber 
needed to produce temperature inversions in the atmo- 
spheres of HD 209458b and XO-lb is currently unknown, 
there appears to be an intriguing connection between the 
equilibrium temperature of the planet in question and the 
presence or absence of an inversion. HD 209458b, which 
has an inversion, receives two times more radiation per 
unit area from its primary than either HD 189733b or 
TrES-1, which do not have inversions. XO-lb is an ex- 
ception to this rule, as it apparently possesses a temper- 
ature inversion despite levels of irradiation comparable 



2 



to those of HD 189733b and TrES-1. Leaving XO-lb 
aside for the moment, it is possible that the increased ir- 
radiation experienced by HD 209458b relative to these 
other planets might have crossed a threshold beyond 
which non-equilibrium compounds form through photol- 
ysis (Burrows et al. 2008), and that these, in turn, pro- 
vide the additional opacity at altitude needed to produce 
the inversion. Similarly, it has been suggested (Hubeny 
et al. 2003; Burrows et al. 2007, 2008; Fortney et al. 2008) 
that the presence of gas-phase TiO or VO in the atmo- 
sphere might also provide the necessary opacity. These 
compounds are predicted to have condensed out of the 
relatively cool atmospheres of HD 189733b and TrES- 
1, whereas they might still remain in gas phase in the 
hottest regions on the day side of HD 209458b. 

We present below observations that allow us to test 
these models by examining the emission spectrum of a 
planet with even higher levels of irradiation than HD 
209458b. TrES-4 orbits a 1.22 Mq star with a period of 
only 3.55 days (Mandushev et al. 2007), and as a result 
it receives twice as much radiation as HD 209458b, and 
four times as much as HD 189733b, TrES-1, and XO-lb. 
Although the TrES-4 primary is more massive than the 
stars in these four systems it is also more distant, and has 
an apparent V-band magnitude of only 11.59 and a K- 
band magnitude of 10.33. The planet TrES-4 is predicted 
to have an equilibrium temperature of approximately 
1700 K, which would make it one of the hottest known 
transiting exoplauets. It also has an unusually large ra- 
dius, even larger than the anomalous radii of HD 209458b 
(Knutson et al. 2007a), WASP-lb (Charbonneau et al. 
2007; Shporer et al. 2007), and TrES-2 (O'Donovan et 
al. 2006), which is difficult to explain without invoking 
additional heating processes such as ongoing tidal circu- 
larization (Bodenheimer et al. 2001, 2003; Baraffe et al. 
2003; Liu et al. 2008). Although the time scale for orbital 
circularization is short compared to the age of the TrES-4 
system, it is possible that interactions with an unknown 
second planet might be pumping TrES-4's eccentricity. 
The current radial velocity data set (Mandushev et al. 
2007), which consists of four points, appears to be con- 
sistent with a circular orbit, although Mandushev et al. 
(2007) do not place an upper limit on the value of the 
eccentricity. Liu et al. (2008) estimate that an eccentric- 
ity of approximately 0.04 could explain TrES-4's inflated 
radius; this would still be consistent with the radial ve- 
locity data. By measuring the timing of the secondary 
eclipse of TrES-4 we will be able to constrain directly 
the planet's orbital eccentricity, either confirming or rul- 
ing out ongoing circularization as the explanation for the 
planet's inflated radius. 

2. OBSERVATIONS 

We observed three secondary eclipses of TrES-4 over a 
period of two weeks in Oct. 2007 using the Spitzer Space 
Telescope (Werner et al. 2004), obtaining data at 3.6, 4.5, 
5.8, 8.0, and 16 /im. We used the Infrared Spectrograph 
(IRS; Houck et al. 2004) on UT 2007 Oct. 8 to observe 
an eclipse in the 16 /im peak-up imaging mode, ac;quiring 
a total of 736 images spanning 7.9 hours with an integra- 
tion time of 30 s for each image. Next we observed a 
secondary eclipse on UT 2007 Oct. 19 using the Infrared 
Array Camera (IRAC) (Fazio et al. 2004), obtaining data 
simultaneously at 3.6 and 5.8 /im. Because this star is 



significantly dimmer than HD 189733 or HD 209458, we 
were able to observe in full array mode using the same 
10.4 s integration time in both channels while still re- 
maining well below saturation in the 3.6 fim channel, 
acquiring a total of 2164 images in each channel over 7.8 
hours. We observed a third secondary eclipse on UT 2007 
Oct. 22 at 4.5 and 8.0 /im, again using 10.4 s integration 
times in both channels and spanning 7.8 hours for a total 
of 2164 images in each channel. 

Because the two shortest wavelength IRAC channels 
(3.6 and 4.5 ^m) use InSb detectors and the three longer 
wavelength channels (5.8, 8.0, and 16 /im) use Si;As de- 
tectors, there are fundamental differences between the 
properties of the data taken with these two types of de- 
tectors. We describe our analysis for each type of detec- 
tor separately below. 

2.1. 3.6 and 4-5 fJ^m Observations (InSb Detector) 

The background in the 3.6 and 4.5 /xm IRAC band- 
passes is extremely low relative to the flux from TrES-4, 
and contributes only 0.25% and 0.30%, respectively, of 
the total flux in an aperture with a 3-pixel radius cen- 
tered on the position of the star. As a result, we found 
that we obtained optimal results using aperture photom- 
etry with a radius of 3.0 pixels to estimate the flux from 
the star in these two channels. We allow the position of 
our aperture to shift with the position of the star in each 
image, and see no evidence for trends correlated with the 
shifting position of the star on the array in the resulting 
time scries for aperture radii larger than 2 pixels. Our 
choice of a 3.0-pixel radius aperture minimizes the like- 
lihood that transient hot pixels will be included within 
that aperture. We determine the position of the star in 
each image as the position-weighted sum of the flux in 
a 5 X 5 pixel box centered on the approximate position 
of the star. We estimate the background in each image 
from an annulus with an inner radius of 12 pixels and an 
outer radius of 50 pixels centered on the position of the 
star. We calculate the JD value for each image as the 
time at mid-exposure and apply a correction to convert 
these JD values to the appropriate HJD, taking into ac- 
count Spitzer's orbital position at each point during the 
observations. 

Fluxes measured at these two wavelengths show a 
strong correlation with the changing position of the star 
on the array, at a level comparable to the depth of 
the secondary eclipse. This effect is due to a well- 
documented intra-pixel sensitivity (Reach et al. 2005; 
Charbonneau et al. 2005, 2008; Morales-Calderon et al. 
2006; Knutson et al. 2008), and can be removed by fitting 
the data with a linear function of x and y position: 

f = fo * {ci + C2{x - xo) + C3{y - yo j) (1) 

where fo is the original flux from the star, / is the 
measured flux, x and y denote the location of the flux- 
weighted centroid of the star on the array, a;o a-nd yo are 
the coordinates of the center of the pixel containing the 
peak of the star's point spread function, and ci — C3 are 
free parameters in the fit. In the 3.6 /xm channel Xq and 
yo had values of [89.5,157.5], and in the 4.5 /itm chan- 
nel they had values of [84.5,156.5]. In contrast to our 
previous observations of HD 189733 and HD 209458 in 
these channels (Knutson et al. 2008; Charbonneau et al. 
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Fig. 1. — Secondary eclipse of TrES-4 observed on UT 2007 
Oct. 19 at 3.6 and 5.8 /im and on UT 2007 Oct. 22 at 4.5 and 
8.0 /im. Data are binned in 6.6 minute intervals and normalized 
to one, then offset by a constant for the purposes of this plot. The 
overplotted curves show the best-fit corrections for detector effects 
(see §2.1 and §2.2). 
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Fig. 2.— Secondary eclipse of TrES-4 observed on UT 2007 
Oct. 19 at 3.6 and 5.8 /im and on UT 2007 Oct. 22 at 4.5 and 
8.0 /jja, with best-fit eclipse curves overplotted. Data have been 
normalized to remove detector effects (see discussion in §2.1 and 
§2.2), and binned in 6.6 minute intervals, then offset by a constant 
for the purposes of this plot. 



2008), we find that adding quadratic terms to this equa- 
tion docs not improve the fit. This is because the total 
drift in x and y positions during the TrES-4 observations 
was a factor of five smaller than for our observations of 
HD 209458 and HD 189733 (see Knutson et al. 2008, for 
a full explanation of the pointing drifts introduced by 
cycling through the subarrays). We find the position of 
the star on the array varied by 0.07 pixels in x and 0.21 
pixels in y during our 3.6 fim observations. During our 
4.5 /um observations the position of the star varied by 
0.12 pixels in x and 0.14 pixels in y. Because the drift 
in X is relatively small during the 3.6 /im observations, 
we find that we obtain the same results if we fit the data 
with a linear function of y position only, removing the 
linear function of x from Eq. 1. We use this simpler 
fit for the 3.6 /zm data in order to reduce the degrees 
of freedom in our fit. In the 4.5 fim channel we obtain 
optimal results using linear functions of both x and y, as 
described in Eq. 1. In both bandpasses the value for 
the fits is not improved by the addition of higher-order 
terms in x and y, and the addition of these higher-order 
terms did not significantly alter the best-fit values for the 
best-fit eclipse times and depths. 

After correcting for the intrapixel sensitivity, a linear 
trend is still visible in both channels with a slope of 
-|-0.030±0.004% per hour at 3.6 /xm and -0.020±0.003% 
per hour at 4.5 /im. There are several possible sources 
of such a linear trend, including (a) variability caused by 
spot activity on the star, (b) the planet's phase curve. 



or (c) a previously uncharacterized instrumental effect 
related to the detector or telescope. The total change 
in fiux over the 8 hours spanned by these observations is 
larger than the secondary eclipse depth in either channel, 
which rules out (b) as an explanation. TrES-4 is a late 
F star and thus should not have significant spot activity; 
the star has a rotation period > 9.3 days (Mandushev et 
al. 2007) and our IRAC observations in these two chan- 
nels were separated by only 3 days, making (b) unlikely 
but not impossible. As a test we perform the same anal- 
ysis on a second bright star in the array with 47% of the 
flux of TrES-4. This star shows the same positive trend 
at 3.6 /xm and negative trend at 4.5 /im after correcting 
for the intrapixel sensitivity. A preliminary analysis of 
similar observations of TrES-2 (F. T. O'Donovan, private 
commun.), which has a comparable brightness to TrES-4, 
also shows the same positive linear trend at 3.6 /im and 
negative linear trend at 4.5 /xm. This argues strongly for 
(c) as the explanation. 

We correct for this previously undocumented effect by 
fitting the data in both channels with a linear function of 
time. This term is fitted simultaneously with the transit 
curve and the correction for the intrapixel sensitivity so 
that we can accurately characterize the additional uncer- 
tainty in the depth and timing of the eclipse introduced 
by these corrections. This means that at 3.6 /xm we are 
fitting for five parameters, including a constant term, a 
linear function of y position, a linear function of time, the 
eclipse depth, and the eclipse time. At 4.5 /xm we fit for 
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the five parameters listed above plus a linear function of 
X position. We also trim the first hour of data from the 
3.6 /um time series for reasons described at the end of this 
section. We fit the data using a Markov Chain Monte 
Carlo method (Ford 2005; Winn et al. 2007) with 10'^ 
steps, where we set the uncertainty on individual points 
equal to the standard deviation of the out-of-transit data 
after correction for the various detector effects. Before 
beginning our fits wc do an intial trim to remove tran- 
sient hot pixels that fell within our aperture, fitting the 
data with either a linear (at 3.6 ^m) or a quadratic (at 
4.5 /im) function of time to remove large trends and then 
discarding outliers lower than 3.5 a or higher than 3.0 tr. 
This corresponds to an effective range of 0.992 — 1.007% 
at 3.6 /um, and a range of 0.988 — 1.010% at 4.5 fim. 
We chose a tighter upper limit for this step because hot 
pixels tend to have values that are too high, rather than 
too low, and this limit provides more effective filtering 
for these pixels. 

Next we carry out the Markov chain fit on the trimmed 
data. We allow both the depth and timing of the sec- 
ondary eclipse to vary independently for the eclipses 
at each of the two observed wavelengths, and take the 
other parameters for the system (planetary and stel- 
lar radii, orbital period, etc.) from Mandushev et al. 
(2007). We calculate our eclipse curve using the equa- 
tions from Mandel & Agol (2002) for the case with no 
limb-darkening. During each step of the chain we ex- 
clude outliers greater than either 3ct (for the 3.6 fim. fit) 
or 4cr (for the 4.5 /zm fit), as determined using the resid- 
uals from the model light curve, from our evaluation of 
the function. We rescale the value of the function 
to account for the fact that we are varying the number 
of pixels included in the fit. Because the correction for 
the intra-pixel sensitivity is larger at 3.6 /im, we find the 
solution in this bandpass is particularly sensitive to the 
presence of points with hot pixels in the time series, and 
we cannot obtain consistent results for the eclipse depth 
over a range of apertures unless we use a tight 3ct limit 
in our fits. With this limit, we obtain consistent results 
for apertures ranging from 3 — 5 pixels. The fit in the 
4.5 /im bandpass is less sensitive to these hot pixels, and 
we achieve consistent results for apertures ranging from 
3 — 5 pixels using a Aa limit. 

After running the chain, we search for the point in 
the chain where the value first falls below the me- 
dian of all the values in the chain (i.e. where the 
code had first foimd the best-fit solution), and discard 
all the steps up to that point. We take the median of the 
remaining distribution as our best-fit parameter, with 
errors calculated as the symmetric range about the me- 
dian containing 68% of the points in the distribution. 
The distribution of values was very close to symmetric 
in all cases, and there were no strong correlations be- 
tween variables. Interestingly, we find that the slope of 
the linear function of time in the 3.6 /im fit increases with 
aperture size; for an aperture with a radius of 4 pixels it 
has increased its slope by 50% relative to that for a 2.5 
pixel aperture. Similarly we find that the slope is steeper 
in the background-subtracted time series for the fainter 
comparison star when compared to TrES-4 photometry 
with the same aperture size. Moreover, the trend in the 
comparison star photometry appears to be more asymp- 
totic than linear in nature, with a steep rise of up to 1% 
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Fig. 3. — Median background flux in the 3.6 array over the 8 
hours spanned by our observations. Note the asymptotic behavior 
of the measured flux values at early times. This same asymptotic 
behavior is also visible in the time series for faint stars visible in 
the 3.6 /im images, albeit with a much smaller (1% instead of 20%) 
amplitude. 

in the first hour and a half of observations and a more 
gradual slope over the rest of the time series. This would 
seem to hint at an illumination-dependent effect, similar 
to the detector ramp described in §2.2 for the 5.8 and 

8.0 /im arrays. 

If this effect behaves similarly to the detector ramp, 
we should see the same trend with a larger amplitude 
in the background flux, which is estimated using pixels 
with lower ilhimination levels than those of either star. 
We find the median background flux in the 3.6 /im band- 
pass increases by 20%i in the first hour of observations, 
and then has a smaller downward trend over the rest 
of the observations (see Figure 3). Although the back- 
ground flux does not have the positive linear slope ob- 
served in the background-subtracted fluxes for TrES-4 
and the comparison star, the increased amplitude (20% 
vs 1%) of the asymptotic rise in the background flux and 
the presence of a similar asymptote in the photometry 
for the comparison star, which is half the brightness of 
TrES-4, argues strongly for an illumination-dependent 
effect. This would also provide a reasonable explanation 
as to why this effect was not detected in 3.6 /tm obser- 
vations of HD 189733 and HD 209458, as these stars are 
approximately a factor of 100 brighter than TrES-4 in 
this bandpass. We elect to trim the first hour of data 
from the time scries at 3.6 /im in order to remove the 
region that appears to exhibit asymptotic behavior. The 
background in the 4.5 /im array has no asymptote and 
appears to be flat at a level of 3% or better, therefore we 
continue to use the entire time series in our flts in this 
bandpass. 

Figure 1 shows the final binned data from these fits 
with the best fit normalizations for the detector effects in 
each channel over plotted, and Figure 2 shows the binned 
data once these trends are removed, with best-fit eclipse 
curves ovcrplotted. Best-fit eclipse depths and times are 
given in Table 1. 

2.2. 5.8 and 8.0 /tm Observations (Si:As Detector) 
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At longer wavelengths the flux from the star is smaller 
and the zodiacal background is larger; we find that the 
background contributes 11% and 12% of the total flux 
in a 2.5-pixel aperture at 5.8 and 8.0 ^m, respectively. 
Because the background is higher in these two chan- 
nels (the median background flux is 1.43 MJy/Sr in the 
5.8 fxm bandpass and 0.71 MJy/Sr in the 8.0 fxm band- 
pass), we used a psf fit to derive the time series in both 
bandpasses and compared the results to those from aper- 
ture photometry. We perform aperture photometry on 
the images in both bandpasses using a radius of 2.5 pix- 
els as described in §2.1. For our psf fits we first calculate 
the background in each image by iterativcly trimming 
3(7 outliers and fitting a Gaussian to the central region 
of a histogram of the remaining pixels. We then sub- 
tract the background from the image and fit the remain- 
ing flux distribution with in-flight point spread functions 
generated from calibration test data^. We use a circu- 
lar region (rounded to the nearest integer pixel) with a 
radius of 3.5 pixels centered on the position of the star 
for our psf fits, and we use the error arrays generated by 
the standard Spitzer pipeline to determine the relative 
weighting for individual pixels. Increasing or decreasing 
the radius of this region by one pixel increased the scat- 
ter in the resulting time series. We also rescale the x and 
y coordinates of our interpolated model psf by a factor 
of 0.8, which effectively reduces the width of the peak of 
the distribution by 20%. This rescaling reduces the 
value for our fits by a factor of two. To fit the observed 
point spread function, we interpolate our model to 520 
times the resolution of the IRAC array (it is already at 
l/8th of the native pixel scale, and we interpolate by an 
additional factor of 65) and then rebin with the psf cen- 
tered at the desired position, which we allow to vary in 
our fits. This allows us to fit for the x and y position of 
the star to a resolution of 1 /520th of a pixel. The scat- 
ter on a representative 30-minute segment of the final 
fitted positions after removing a linear trend with time 
is ±0.02 pixels, ten times larger than the minimum pixel 
resolution in our fits, so this is a reasonable choice. We 
also fit for a constant scaling factor corresponding to the 
total fiux. 

We flag bad pixels marked by the Spitzer pipeline 
in our subarray and give them zero weight in our fits. 
To find transient hot pixels, we divide the time series 
into sets of 20 images and calculate the median value 
and standard deviation at each individual pixel position 
within each set of 20 images. Wc then step through the 
individual images and mark outliers more than 3(t away 
from the median value for that pixel position as bad pix- 
els in that image. We find that 98% of our images have 
one or fewer bad pixels in the region used for our fits, 
which contains 29 pixels in total. This process signifi- 
cantly reduces the number of large outliers in the final 
time series, although it does not eliminate such outliers 
completely. We found that increasing our threshold for 
bad pixels to 4(t and then lOa outliers produced compa- 
rable results with an increasing number of large outliers 
in the final time series. 

At 5.8 /Ltm we found that the relative scatter in the time 
series from the psf fits was 20% higher than in the time 
series from aperture photometry with a 2.5 pixel radius. 

Available at http://ssc.spitzer.caltech.edu/ira<;/psf.html 



As a result of this increased scatter, which is likely pro- 
duced by discrepancies between the model psf and the 
observed psf, we conclude that aperture photometry is 
still preferable in this channel. We compare the time se- 
ries using apertures ranging from 2 — 3.5 pixels and find 
consistent results in all cases, but with a scatter that 
increases with the radius of the photometric aperture. 

At 8.0 /zm wc found that psf photometry produced a 
time series with an out-of-transit RMS variation that was 
2% lower than that of the equivalent time series using 
aperture photometry with a radius of 2.5 pixels. More 
significantly, the use of a psf fit with bad pixel filtering ef- 
fectively corrected the fiuxes in 30 images with transient 
hot pixels in the photometric aperture, allowing them to 
be used in the final time series. In light of these two 
changes, we conclude that psf photometry is preferable 
in this bandpass. 

After determining the optimal method for estimating 
the fluxes in each bandpass, we must remove any de- 
tector effects in order to determine the best-fit eclipse 
depths. There is no intra-pixel sensitivity in these two 
bandpasses, but there is another well-documented de- 
tector effect (Knutson et al. 2007b, 2008; Charbonneau 
et al. 2008) that causes the effective gain (and thus 
the measured flux) in individual pixels to increase over 
time. This effect has been referred to as the "detector 
ramp" , and has also been observed to occur in the IRS 
16 /im peak-up array, which is made from the same mate- 
rial (Deming et al. 2006). The size of this effect depends 
on the illumination level of the individual pixel. Pixels 
with high (>250 MJy Sr~^ in the 8 /zm channel) will 
converge to a constant value within the first hour of ob- 
servations, whereas lower-illumination pixels will show 
a linear increase in the measured flux over time with a 
slope that varies inversely with the logarithm of the illu- 
mination level. In our observations of TrES-4, this effect 
produces a 1.5% increase in the measured flux from the 
star at 8.0 /xm during the 8 hours spanned by these ob- 
servations (see Figure 1), and a much smaller (0.1%) in- 
crease in the measured flux from the star at 5.6 iim over 
the first hour of observation. In both bandpasses the 
ramp has an asymptotic shape, with a steeper rise in the 
first 30 minutes of observations. We correct for this ef- 
fect by fitting our time series in both bandpasses with 
the following quadratic function of ln{dt): 

f = fo*{ci+ C2ln{dt + 0.02) + cs{ln{dt + 0.02))^) (2) 

where /o is the original fiux from the star, / is the 
measured flux, and dt is the elapsed time in days since 
the start of the observations. In previous observations 
(Knutson et al. 2008; Charbonneau et al. 2008) we 
trimmed the first 30 minutes of data from the time series 
in order to avoid the steepest part of the ramp. This was 
not necessary here, as the ramp is not as steep for fainter 
sources, and the addition of a constant term of 0.02 days 
in Eq. 2 ensures that we are fitting the same function to 
these data as before. It also ensures that wc avoid the 
singularity at d.t = 0. Because TrES-4 is approximately 
100 times fainter than HD 209458 and HD 189733 in this 
bandpass, the pixels at the center of the star's psf receive 
lower levels of illumination and the slope of the ramp 
at early times is correspondingly more gradual. Offset- 
ting the zero point of the curve reduces the slope of the 
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TABLE 1 

Best-Fit Eclipse Depths and Times 



A (/im) 


Eclipse Depth 


nn a 
bright 


Center of Transit (HJD) 


O-C (min.)'' 


3.6 


0.137 ±0.011% 


I960 ± 70 K 


2454392.6137 ± 0.0025 


5.6 ±3.6 


4.5 


0.148 ±0.016% 


1800 ± 90 K 


2454396.1665 ± 0.0032 


4.0 ±4.5 


5.8 


0.261 ± 0.059% 


2210 ± 300 K 


2454392.5942 ± 0.0110 


-22.4 ± 15.8 


8.0 


0.318 ± 0.044% 


2290 ± 220 K 


2454396.1696 ±0.0045 


8.4 ±6.5 


16.0 


< 0.37%"^ 


< 2350 K° 







^We calculate the planet's brightness temperature in each bandpass as follows: first we set 
the planet's emission spectrum equal to a Planck function with the given temperature, then 
we take the flux-weighted average of the planet-star flux ratio over the bandpass in question, 
solving for the temperature required to match the observed eclipse depth in that bandpass. 
We use a 6200 K Kurucz atmosphere model for the stellar spectrum; this is the same stellar 
spectrum used to calculate the planet-star flux ratios plotted in Fig. 6. 

^Observed minus calculated transit times, where the expected transit times are calculated 
using the ephemeris from Mandushev et al. (2007) and assuming zero eccentricity. 

^ These are the 2a upper limits on the eclipse depth and brightness temperature in this 
channel. 



function at early times and provides a better fit to the 
observed behavior. We note that our specific choice of 
constant has a ncgHgiblc effect on the final vahics for the 
echpsc depths in these two channels; repeating these fits 
with a constant of O.Of or 0.03 changes the best-fit eclipse 
depth at 5.8 /xm by O.OScr and the best-fit eclipse depth 
at 8.0 fim by 0.06fT. 

We fit both Eq. 2 and the transit curve to the data si- 
multaneously using a Markov Chain Monte Carlo method 
as described in §2.1. As before, the distribution of val- 
ues was very close to symmetric in all cases, and there 
were no strong correlations between the variables. Best- 
fit eclipse depths and times from these fits are given in 
Table f , and the time series before and after correcting 
for detector effects are shown in Figures 1 and 2, respec- 
tively. As a check we repeated these fits without the 
quadratic term in Eq. 2, and found that the value of the 
function for our best-fit solution increased by 0.3 and 
3.3 at 5.8 and 8.0 /im, respectively. Our choice of a linear 
or quadratic function has a negligible effect on the best- 
fit eclipse depth in the 5.8 /xm bandpass, therefore we 
opt for a more consistent approach and retain the same 
quadratic function for our fits in both bandpasses. 

2.3. 16 fxm Observations (Si: As Detector) 

The median background in the 16 fxm IRS peak-up ar- 
ray during our observations is 1.85 MJy/Sr, more than 

twice the level of the background in the 8.0 /im band- 
pass, while the star is correspondingly fainter. In this 
bandpass we find that the background contributes 71% 
of the total flux in a 2.0 pixel aperture, and it is the 
primary source of noise in our final time series. We es- 
timate the flux from the star in each image using both 
aperture photometry with a radius of 2 pixels (which 
minimizes the RMS scatter in the resulting time series 
for aperture photometry) and a psf fit with a radius of 
2.5 pixels, and then compare the results. We estimate 
the sky background using the central 33 x 44 pixel region 
of the array, iteratively trimming 3a" outliers and fitting 
a Gaussian to the central region of a histogram of the 
remaining pixels. Excluding the region containing TrES- 
4 from this histogram changed the median background 
value by only 0.02%, a negligible amount. For our psf fit 
we calculate the median background for each image and 
subtract that value, then fit the remaining flux distri- 



bution with a model point spread function derived from 
observations of HD 42525 in this bandpass* , where we al- 
low the X and y position of the model to vary freely, along 
with a constant scaling factor corresponding to the total 
flux. We use the error arrays generated by the standard 
Spitzer pipeline to determine the relative weighting for 
individual pixels. As described in §2.2, we rescale the x 
and y coordinates of our interpolated model psf by a fac- 
tor of 0.9, which effectively reduces the width of the peak 
of the distribution by 10%, in order to provide a better 
match for the width of the star's psf in our images. We 
find that increasing or decreasing the 2.5 pixel radius of 
the region used for our psf fits increases the scatter in the 
resulting time scries. In light of the higher background 
fluxes in this channel, it is not surprising that we obtain 
better results with a smaller region than the 3.5-pixel re- 
gion used in our 8.0 /xm fits. We find that the relative 
RMS scatter in the final time series is 17% lower for psf 
fits than for aperture photometry, and we use the psf fits 
for all of our subsequent analysis. 

We plot the time series for the psf fits in the upper 
panel of Figure 4; as this figure illustrates, the point- 
to-point scatter in the time series is high enough that 
we are unable to detect the secondary eclipse. We esti- 
mate the RMS variation in the time series by fitting the 
data outside of the predicted time of eclipse with a linear 
function of time, dividing by this function and then iter- 
atively trimming 3(t outliers. The standard deviation of 
the remaining points is 3.3%. We compare this to a sim- 
ple estimate! of the predicted noise in the 2 pixel aperture 
that we used for our aperture photometry, including the 
read noise, dark current, and photon noise from the star 
and background, as described in §7.2.3.4.4. of the Spitzer 
Observer's Manual: 



N = 



{isky ~\~ idark ~\~ istar) * tint "I" 2^77 



{RNf 



reads 



(3) 



where N is the total noise in electrons, isky and istar 
are the electrons generated per second in our two-pixel 
aperture from the sky background and the stellar flux, 
respectively, idark is the dark current contribution, tint 

^ Available at http://ssc.spitzer.caltech.edu/irs/puipsf/. 
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Fig. 4. — Data from our 16 ^in observations spanning a predicted 
time of secondary eclipse on UT 2007 Oct. 8. The upper panel 
shows the time series with the best-fit correction for the detector 
ramp overplottcd (sec discussion in §2.3). The lower panel shows 
the data after it has been normalized to remove this ramp, with the 
best-fit eclipse light curve overplotted. This light curve increases 
slightly in flux during the eclipse event, although this increase is 
not statistically significant. Data in both panels has been binned 
in 6.5 minute intervals. 



is the integration time, RN is the read noise in e~ 
pixel"^, and Nreads is the number of reads. Our images 
have a total integration time of 31.46 s, 16 reads per im- 
age, and we use a value of 30 e^ pixel^^ for the read 
noise. The median background flux during our observa- 
tions is 1873 e~ s~^ in our two pixel aperture, the dark 
current is < 126 e~ s~^, and the star contributes a me- 
dian flux of 777 c~ s~^ in this aperture. Evaluating Eq. 
3 with these numbers and assuming the dark current is 
negligible, we find a total noise contribution of 290 e~ in 
each 31.46 s exposure, which would translate to a relative 
RMS of 1.2%. 

Compared to this simple estimate, we find that the rel- 
ative RMS variation in the time series of 3.3% is a factor 
of 2.8 higher than expected. This is surprising, but an ex- 
amination of the scatter in the fluxes for individual back- 
ground pixels across the images indicates that this back- 
ground has a noise level three times higher than predicted 
based on photon and read noise alone. This increased 
background noise could easily explain the increased scat- 
ter in the final time series. Deming et al. (2006) ob- 
served a secondary eclipse of HD 189733b at 16 fjLm and 
found that they were able to obtain a RMS variation 1.7 
times higher than the predicted photon noise, but this 
star is approximately 100 times brighter than TrES-4 in 
the infrared and the relative noise contribution from the 
sky background is correspondingly small in these obser- 
vations. Although the RMS variation in the backgrotmd 
fluxes for our TrES-4 observation is inconsistent with the 
predicted photon and read noise alone, it is consistent 
with the uncertainties produced by the Spitzer pipeline. 
The Spitzer error estimates include additional uncertain- 
ties from dark current subtraction, droop correction, flat 
flelding, and other other steps in the standard Spitzer 
pipeline, therefore we speculate that it is one (or more) 



0) 

Q 



o 



0.02- 



0.01 - 



0.00 




-0.5 0.0 
Depth of Eclipse ( 



Fig. 5. — Probability distribution for the eclipse depth from a 
Markov Chain Monte Carlo fit to the IRS data. Negative values 
for the eclipse depth correspond to an increase in flux during the 
eclipse event, and positive values correspond to a decrease in flux 
(as would normally be expected for an eclipse). The long-dashed 
and short-dashed lines indicate the 2a and Scr limits on the eclipse 
depth, respectively, which were calculated by integrating over this 
distribution. 

of these steps that is the dominant source of noise in our 
data. 

We place an upper limit on the eclipse depth in 
this bandpass by fitting the data simultaneously with 
a quadratic function of ln{dt) as described in Eq. 2 and 
an eclipse function where we have fixed the timing to the 
predicted value and allow the depth to vary over both 
positive and negative values. We fit the data using a 
Markov chain Monte Carlo method as described in §2.1, 
with uncertainties for individual points set equal to the 
RMS variation in the out-of-transit data after the de- 
tector ramp has been removed. Because the star is not 
much brighter than the background, pixels at the cen- 
ter of the star's psf exhibit a similar ramp to that of 
the background pixels, and subtracting this background 
removes the majority of the detector ramp described in 
§2.2. There is still a small ramp remaining in the final 
time series (see upper panel of Figure 4), and we find 
that fitting this ramp with Eq. 2 reduces the value 
of the best-fit solution by 0.5 when compared to a sim- 
ple linear function of time. Although this is not a large 
improvement, we note that the quadratic fit is a more 
accurate description of the detector ramp, and using it 
in our fits here allows to determine the upper limit on the 
eclipse depth in a manner consistent with our previous 
analysis in the 5.8 and 8.0 ^m bandpasses. 

Our best-fit eclipse depth from these fits is —0.14 ± 
0.30%, indicating that the measured flux increased by 
a statistically insignificant amount during the eclipse 
event. Integrating over the probability distribution for 
this parameter from the Markov fit (see Fig. 5), we find 
that the 2cr upper limit on the eclipse depth is 0.37%, 
and the Sa upper limit is 0.69%. We adopt the 2a limit 
in the discussion below. 

3. DISCUSSION 

We find that the RMS variation in our final time series 
is 0.96, 1.06, 1.24, and 1.35 times the predicted photon 
noise from the star and background flux at 3.6, 4.5, 5.8, 
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Wavelength (microns) 

Fig. 6. — Day-side planct-star flux ratios for TrES-4 as determined from measurements of the secondary eclipse depth in the four IRAC 
bandpasses (red circles). The horizontal bar and arrow at 16 fim show the 2a upper limit on the eclipse depth in this bandpass. The 
black line corresponds to the default model (no temperature inversion) with a redistribution parameter P„ = 0.1, which describes the 
case where 10% of the incident energy is redistributed to the night side. The purple, blue, and green lines correspond to models with an 
additional optical absorber at high altitudes (parameterized as Kextra), which produces a thermal inversion around pressures of 0.001 bar 
(Burrows et al. 2007, 2008). Squares show the values for these models after integrating over the Spitzer bandpasses. The high planet-star 
flux ratios at 5.8 and 8.0 ^m argue strongly for the presence of an inversion, as models with K^xtra > cw? /g provide the best match at 
these wavelengths. The eclipse depths in the 3.6 and 4.5 /im bandpasses are best matched by a model with relatively efficient day-night 
circulation and modest additional opacity (purple model, K^xtra = 0.1 cm-^/g and P„ = 0.3). 



and 8.0 /im, respectively. In the 16 /im bandpass the 
noise is a factor of 2.8 higher than the predicted contri- 
bution from the photon noise of the star and background 
and the detector read noise, but it is consistent with the 
estimated uncertainties from the Spitzer pipehne. As a 
resuh of this increased uncertainty, we are unable to de- 
tect the secondary eclipse in the 16 /im bandpass, and 
instead place an upper limit on the eclipse depth at this 
wavelength. 

We determine the best-fit eclipse times for the two 
secondary eclipses observed using IRAC by taking the 
weighted average of the best-fit eclipse times in each 
bandpass. Using this method, we find that the eclipse 
observed in the 3.6 and 5.8 /im bandpasses occurred 
4.2 ± 3.5 minutes later than the predicted time based 
on the ephemeris from Mandushev et al. (2007), where 
we have neglected the light travel time in the TrES-4 sys- 
tem (on the order of 30 s; Loeb 2005) and assumed the 
secondary eclipse will occur exactly half an orbit after 
the transit. We repeat this analysis for the eclipse ob- 
served in the 4.5 and 8.0 /im bandpasses, and find that it 
occurred 5.5 ± 3.7 minutes later than predicted. If we as- 
sume that the planet's orbit remained the same over the 
3.5-day period spanned by our observations (ie no per- 
turbations that would change the orbital semi-major axis 
or eccentricity during this time), we can combine obser- 
vations in all four of the IRAC bandpasses to get a single 
estimate for the best-fit eclipse time, which we find occurs 
4.8 ±2.6 minutes later than predicted. However, there is 
an additional ±5.0 minute uncertainty in the predicted 
transit time from Mandushev et al. (2007). Including 
this uncertainty we find the two averaged eclipses occur 



4.8±5.6 minutes later than predicted, which is consistent 
with zero offset. 

Our estimate for the best-fit timing offset translates 
to a constraint on the orbital eccentricity e and the ar- 
gument of pericenter w of ecos(w) = 0.0015 ± 0.0017; 
the 3(7 upper limit on this value is |ecos(w)| < 0.0058, 
where we have calculated the limit by integrating over 
a Gaussian distribution with limits of integration that 
are symmetric around zero. We selected these limits be- 
cause we are interested in constraining the magnitude of 
e rather than the sign of the cos(aj) term. This upper 
limit means that unless the longitude of periastron lj is 
close to 90° or 270°, we can rule out tidal heating from 
ongoing orbital circularization (Bodenheimer et al. 2001; 
Liu et al. 2008) as an explanation for TrES-4's inflated 
radius. Liu et al. (2008) estimate that TrES-4 would need 
to have an orbital eccentricity of approximately 0.04 to 
provide the required energy; this would require 1^1 to be 
less than 9° away from the two angles listed above to be 
consistent with our upper limit. Our conclusion is also 
consistent with fits to the four radial velocity points from 
Mandushev et al. (2007), although these points provide 
a relatively weak constraint on the eccentricity. 

Next we compare the secondary eclipse depths in the 
four IRAC bandpasses to the predictions from atmo- 
sphere models for this planet (see Fig. 6 and 7 below). 
In order to fit the IRAC data and the 16-micron up- 
per limit, we employed the same formalism described in 
Burrows et al. (2007, 2008). Using the planet-star radius 
ratio of 0.09903 ± 0.00088 from Mandushev et al. (2007) 
and a Kurucz atmosphere model (Kurucz 1979, 1994, 
2005) with an effective temperature of 6200 ± 75 K for 
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Fig. 7. — Day-side pressure-temperature profiles for the four 
models plotted in Fig. 6. Increasing the opacity of the optical 
absorber (parameterized as Kextra) increases the temperatures at 
pressures less than 0.01 — 0.001 bars, while increasing the fraction of 
heat redistributed to the night side (parameterized as P„ ) decreases 
the temperatures around pressures of 0.01-0.1 bars. This is because 
Burrows et al. (2008) parameterize the effects of energy transport 
to TrES-4's night side in these ID models by adding a heat sink at a 
pressure of 0.1 bars, which causes a drop in temperature as the day- 
night circulation is turned up. See Burrows et al. (2007, 2008) for 
a full description of these parameterizations and the corresponding 
models. 

the stellar spectrum (Sozzetti et al. 2008), we calculated 
the emergent spectrum at secondary eclipse for a pair of 
free parameters, P„ and Kextra- Pn is the dimensionless 
redistribution parameter that, in approximate fashion, 
accounts for the cooling of the dayside and the warming 
of the nightside by zonal winds near an optical depth of 
order unity. It is a measure of the efficiency of heat re- 
distribution by super-rotational hydrodynamic flows. As 
the value of P„ is increased, the day side becomes cooler 
and the emergent planetary flux at superior conjunction 
becomes correspondingly small. Kextra is the absorptive 
opacity in the optical at altitude (here in cm^/g) used 
to create a temperature inversion. The origin of such an 
absorber is unknown, but might be due to exotic species 
generated by the severe insolation regime or to TiO/VO 
at altitude (Hubeny et al. 2003; Burrows et al. 2007, 2008; 
Fortney et al. 2008). Concerning the latter, the coldtrap 
effect should deplete the upper atmosphere of such di- 
atomics, but atmospheric circulation and/or turbulence 
could alter this picture by advecting these compounds 
up from the lower atmosphere or removing them on the 
cooler night side of the planet. 

Fig. 6 and 7 show four models with varying values 
for Pn and Kextra- The standard non- inverted model 
{Kextra — cm^/g) is clcarly inconsistent with the ob- 
served fluxes from TrES-4 at wavelengths longer than 
4 /im. It is possible to match the observed 3.6 //m flux 
with this model by reducing the relative fraction of the 
incident energy that is redistributed to the planet's night 
side, thus increasing the day-side temperature and cor- 
responding fluxes, but even this change is insufflcient at 
longer wavelengths. In contrast to this model, all three 
models with a thermal inversion {Kextra > cm^/g) pro- 
vide an improved match to the 5.6 and 8.0 fim fluxes. 



The ratio of the 3.6 fiui and 4.5 /im fluxes, another 
measure of inversion (Burrows et al. 2007), is also less 
than one, much lower than for models without inversions. 
The best overall fit is obtained by setting P„ = 0.3 and 
Kextra = 0.1 cm^/g. Corresponding to a case with rel- 
atively efficient day-night circulation and modest addi- 
tional opacity. 

We note that our particular choice of planet-star ra- 
dius radio and stellar effective temperature may affect 
the predictions of these models. As a test we re-run 
our full radiative transfer codes for the P„ = 0.3 and 
Kextra — 0.1 cm^/g modcl using stellar atmosphere mod- 
els with temperatures of 6100 K and 6300 K, and mea- 
sure the resulting change in the predicted eclipse depths 
in the four IRAC channels. We find that decreasing the 
stellar effective temperature by 100 K shifts the pre- 
dicted eclipse depths in the [3.6,4.5,5.8,8.0] fim band- 
passes by [-1.8%,-1.3%,-l.l%,-0.9%], while increas- 
ing the temperature by 100 K results in changes of 
[+1. 7%, -Hi. 4%, -hi. 1%, +0.9%] in these same bandpasses. 
These changes are negligible relative to the la uncertain- 
ties of [8.0%, 10.8%,22.6%, 13.8%] in the measured eclipse 
depths. From this test we conclude that increasing or 
decreasing the effective temperature of the star by an 
amount comparable to the formal uncertainties in this 
quantity cannot produce a non-inverted model that is 
consistent with the measured eclipse depths at longer 
wavelengths, and is unlikely to alter our conclusion that 
the Pn — 0.3 and Kextra — 0.1 cm^/g model provides the 
best fit to the data. The predicted eclipse depths scale 
linearly with the planet-star area ratio, but squaring the 
planet-star radius ratio from Mandushev et al. (2007) 
results in a value of 0.98 ± 0.02% for this quantity; this 
is a factor of 100 smaller than the typical uncertainty 
contributed by the stellar effective temperature. 

It is interesting to note that the same values of Kextra = 
0.1 cm^/g and P„ = 0.3 also provide the best fit to the 
observed broadband emission spectra for HD 209458b 
(Burrows et al. 2007, 2008). This would seem to im- 
ply that the atmospheric circulation and relative abun- 
dances of the species responsible for creating the inver- 
sions in the upper atmospheres of both planets are simi- 
lar, despite the higher temperatures and increased levels 
of irradiation experienced by TrES-4. This is probably 
an oversimplification of the problem, however, as there 
are likely substantial thermal and chemical gradients be- 
tween the substellar point and the day-night termina- 
tor on both planets, and our observations constrain only 
the hemisphere-averaged properties of the day-side at- 
mosphere. 

4. CONCLUSIONS 

Our observations of TrES-4 at 3.6, 4.5, 5.8, 8.0, and 
16.0 /im reveal that this planet has a thermal inversion 
similar to the one observed for HD 209458b (Knutson 
et al. 2008; Burrows et al. 2007). The presence of an 
inversion in the atmosphere of TrES-4 provides support 
for the idea that planets with higher levels of irradiation 
are more likely to have thermal inversions, although it 
does not distinguish between competing theories for the 
nature of the optical absorber responsible for the creation 
of the inversions. 

If we are to fully understand the mysterious origin 
of these temperature inversions, it will require a much 
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larger sample than the seven systems with published 
Spitzer observations. There are currently 40 known tran- 
siting planetary systems with published coordinates, of 
which 33 are bright enough to observe with Spitzer. A 
sample spanning a range of stellar metallicities, levels of 
irradiation, surface gravities, and orbital periods might 
reveal important connections between the presence of a 
temperature inversion and other properties of the system. 
Such a comprehensive survey has the potential to provide 
an explanation for XO-lb's thermal inversion (Machalek 
et al. 2008), which is difficult to understand within the 
current irradiation-dependent picture. 

We estimate an upper limit of \ecos{uj)\ < 0.0058 for 
the orbital eccentricity, consistent with the fits to the ra- 
dial velocity data presented by Mandushev et al. (2007). 
This upper limit means that unless the longitude of pe- 
riastron oj is close to 90°or 270°, we can rule out tidal 
heating from ongoing orbital circularization at the level 
required by Liu et al. (2008) in order to explain TrES-4's 
inflated radius. 

Although Spitzer will deplete its supply of cryogen in 
spring 2009, the 3.6 and 4.5 /im channels will continue to 
function at full sensitivity, and a two-year mission using 
these two channels was recently approved by NASA. Ob- 
servations of the secondary eclipse in these two channels 
should be sufficient to distinguish between planets with 
and without temperature inversions, as demonstrated by 



our discussions above. These observations can be com- 
pared to the set of more than a dozen planets for which 
there will be secondary eclipse observations in all four 
IRAC channels by the end of the cryogenic Spitzer mis- 
sion. This basis set of comprehensive observations, when 
combined with a significantly larger survey in the 3.6 and 
4.5 /im channels, should provide the statistical leverage 
needed to resolve the origin and nature of these inver- 
sions. 
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